Viscous coalescence of droplets: a Lattice Boltzmann study 



M. Gross and I. Steinbach 
Interdisciplinary Centre for Advanced Materials Simulation (ICAMS), 
Ruhr-Universitat Bochum, Universitaetsstr. 90a, 44^89 Bochum, Germany 



D. Raabe 

Max-Planck Institut fur Eisenforschung, Max-Planck Str. 



1, 40237 Diisseldorf, Germany 



o 

(N 



(N 

i 

q 

O 

•i-H 

CZ5 

^ ■ 

(N : 
> . 
\o ■ 
m 

in 
o 



F. Varnilfl 

Interdisciplinary Centre for Advanced Materials Simulation (ICAMS), 
Ruhr-Universitat Bochum, Universitaetsstr. 90a, 44789 Bochum, Germany and 
Max-Planck Institut fur Eisenforschung, Max-Planck Str. 1, 40237 Diisseldorf, Germany 

The coalescence of two resting liquid droplets in a saturated vapor phase is investigated by Lattice 
Boltzmann simulations in two and three dimensions. We find that, in the viscous regime, the bridge 
radius obeys a i 1//2 -scaling law in time with the characteristic time scale given by the viscous time. 
Our results differ significantly from the predictions of existing analytical theories of viscous coales- 
cence as well as from experimental observations. While the underlying reason for these deviations 
is presently unknown, a simple scaling argument is given that describes our results well. 



I. INTRODUCTION 

Coalescence of liquid droplets is important in many 
technological and natural phenomena, as, for example, 
coating and sintering processes phase separation of 
emulsions 2], or the formation of rain drops in clouds 
Q . Coalescence is initiated when two droplets come into 
contact and form a liquid bridge, which then starts to 
grow due to surface tension. This growth is typically 
either opposed by viscous dissipation or inertial forces, 
until finally the two droplets have merged to a single 
droplet. 

Assuming that the initial growth of the liquid bridge 
just results from a competition between surface tension 
a and fluid viscosity ry, the characteristic velocity scale 




FIG. 1: Sketch of two coalescing droplets. Ro is the droplet 
radius, b the radius of the connecting bridge and r a the radius 
of curvature of the meniscus. The right image is a magnifica- 
tion of the bridge region. 
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is given by the capillary velocity, u c = cr/rj. Taking the 
relevant length scale to be the size b of the liquid bridge 
connecting the two droplets (see Fig.[T]), a Reynolds num- 
ber can be defined as Re = pu c bjr\ = pab/rj 2 . Due to the 
presence of b, this number will always be small in the 
beginning of the process, regardless of the value of a and 
r\. The domain Re < 1 then defines the viscous regime, 
which can be described by the Stokes equations, in con- 
trast to the inertial regime, Re > 1, where the Eulcr 
equations hold. 

Coalescence in the inertial regime is quite well under- 
stood nowadays P4Io| and it has been established the- 
oretically and experimentally that the evolution of the 
bridge radius in the inertial stage of coalescence follows 
a scaling law of the form b oc {Rqu / pY^t 1 / 2 , where Rq 
is the initial radius of each drop. 

For the viscous regime, however, there remain some 
discrepancies between experiments and theory. The first 
effort to describe viscous coalescence dates back to the 
work of Frenkel [l2[ on viscous sintering of solid parti- 
cles. By considering the balance between reduction of 
surface free energy and viscous dissipation, he derived a 
simple scaling law for the growth of the bridge radius, 
b cx (R a/'q) 1/2 t 1 / 2 , which was confirmed in subsequent 
sintering experiments [l3| . More sophisticated analysis of 
liquid droplet coalescence in the Stokes regime by Hopper 
[Tij and Eggers et al. [15j . aiming to overcome some of 
the simplifications of Frenkel's model, arrived at a char- 
acteristic t am - or ilogi-growth law for the bridge radius 
at early times in the two- and three-dimensional case, re- 
spectively. However, subsequent experiments [HH, BEII 
neither reproduced Frenkel's nor Eggers' results, but in- 
stead reported a linear time dependence, b oc at/rj in the 
viscous regime fill ]. The disagreement between theory 
and experiments was usually attributed to the asymp- 
totic nature of the analytical results, which are expected 
to hold only in the limit of very small bridge radii. In 
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a recent work, Paulsen et al. [17j address these discrep- 
ancies and argue that most experiments have in fact not 
operated in the actual Stokes regime, but instead in a 
different, "incrtially-limited" viscous regime where the 
drops are governed by a balance between surface tension, 
viscous forces and the center-of-mass inertia. 

In this work, we study viscous coalescence of two 
droplets by Lattice Boltzmann computer simulations of a 
single-component, isothermal, non-ideal fluid in two and 
three dimensions. We consider two initially quiescent 
liquid droplets in a surrounding saturated vapor phase 
for various values of surface tension, viscosity and vapor 
density. Interestingly, neither a growth of the bridge ra- 
dius that is linear in time, as reported in experiments 
0, IE S [Hj]> nor a growth as predicted by theories of 
Stokesian coalescence [lj, HU (i.e., b ~ t - 86 in 2D or 
b ~ — tlogt in 3D) is observed by us. Instead, we find 
that the bridge radius is governed by an approximate 
t 1 ' 2 -scaling law, with the characteristic scale factor set 
by the viscous time - similar to Frenkel's original result. 
While we generally find evaporation and condensation 
processes taking place in our system, these are strongly 
suppressed compared to advection and do not seem to 
have a significant impact on the time evolution of the 
bridge radius. Based on these findings, we provide a sim- 
ple scaling argument that describes our results well. 



II. SIMULATION MODEL 
A. Model 

Numerically, we solve the Navier- Stokes equations for 
a non-ideal fluid in both two and three dimensions with 
the Lattice Boltzmann (LB) method, which is a well- 
established tool for fluid dynamical simulations ranging 
from the micro- to the macroscale [2(|. Motivated by 
the growing interest in microfluidic applications, the LB 
method has been widely applied to the study of problems 
encompassing fluids in complex geometries or with mul- 
tiple phases, such as droplet spreading or wetting [2lj |. 
In the present work, we employ an implementation of 
the LB model described in Q, which features faithful 
representation of isothermal two-phase thermo- and hy- 
drodynamics [22l |23|. 

The LB method is based on an evolution equation for 
the distribution function /(r, c) representing the proba- 
bility to find a fluid element at a given position r moving 
with a certain velocity Cj, 

fi(T+Ci,t+l) = fi(v,t)-- [ft(T,t) - fP(r,t)]+f t F (r,t) . 

T 

(1) 

Here, /,(r) = /(r, Cj) stands for the discretized distribu- 
tions living on the computational lattice domain. The 
lattice nodes arc linked by a set of velocity vectors Ci 
(in the present case, we have i = 1,...,9 in 2D and 
i = 1, . . . , 19 in 3D). Furthermore, /° q is a given equi- 
librium distribution and ff is a forcing-term that gives 



rise to a physical body force F in the Navier-Stokes equa- 
tions (see eq. (O below). In the present case, this body 
force mediates the non-ideal fluid effects. For the specific 
forms of f^ q and ff as well as details of the implementa- 
tion we refer to the original publication [8| . The physical 
observables density p and flow velocity u can be obtained 
at any time-step from the lowest-order moments of /, 
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pu = ^/ l c J + iF = ^/rc 4 , 



(2) 



(3) 



As a consequence of the LB dynamics of the present 
model, the density p and velocity u fulfil the Navier- 
Stokes equations of an isothermal, compressible non-ideal 
fluid, which consist of a continuity equation 



d t p + V • (pu) = , 
and a momentum equation, 
dt(pu) + V(puu) = F + ?yV 2 u + (C + [1 



(4) 



2/d\t]) VV-u. 

(5) 

Here, 77 and £ are the shear and bulk viscosity, which 
can be specified as input to the simulation. The body 
force F is responsible for the non-ideal character of the 
fluid and can be derived from a free energy functional. 
The thermodynamics of a non-ideal fluid is governed by 
a Ginzburg-Landau free-energy functional T of the form 



T\p\ = / dr 



UP) 



(6) 



with k being a square-gradient parameter and /o a Lan- 
dau free energy density. Here, we take /o to be a sim- 
ple double- well potential with freely prescribable minima 
around the equilibrium densities pv, Pl [1, H3| 

Mp) = (3(p-pv) 2 (p-Pl) 2 . (7) 

(3 is a free parameter controlling the compressibility of 
the model. From T we can define a generalized chemical 
potential 

dI d p f - kV 2 P (8) 
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and finally the driving force that enters the Navier-Stokes 
eq. © 



F = -pVfi. 



(9) 



In equilibrium, the free energy functional eq. (JB]) admits 
for coexistence of liquid and vapor phases separated by 
a diffuse interface [251 ] . In the case of a flat interface 
located at position x = 0, the density profile is given by 



Pint 

where 



1 1 / x \ 
0) = t;(Pl + Pv) + t;(Pl - /v)tanh (-) , 

2 2 \w/ 



Pl - Pv 

is a measure for the interface width. 



(10) 



(11) 
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B. Setup 

In the 2D case, we study two resting droplets of iden- 
tical radii of R = 5000 lattice units (l.u.), initially sepa- 
rated by 4 l.u., that are placed into a periodic rectangular 
domain of size 20100 x 10200 l.u. For comparison, we also 
performed simulations using droplet radii of Rq = 400 
and box sizes of 1800 x 1000 l.u., but keeping all the other 
simulation parameters the same as in the case of large ra- 
dius. The discussion of the velocity field in sec. MI Al is 
based on simulations of droplets with radii of Rq = 1600. 
Due to computational limitations, in the 3D case simula- 
tions generally have system sizes of 460 x 250 x 250, where 
two spherical droplets each of radius Ro = 100 l.u. and 
initially separated by 4 l.u. are placed. The initial sepa- 
ration of the two droplets has been varied in a few cases, 
showing that this parameter has a negligible influence on 
our results. Also, in a number of cases, simulations with 
larger system sizes but identical droplet radii have been 
performed to ensure that periodic boundary conditions 
have negligible influence on the evolution of the bridge 
radius. The bridge radius b is determined from the (inter- 
polated) position along the horizontal symmetry axis of 
the simulation box that corresponds to a density value of 
{pl+Pv)/2- We follow the time evolution until the bridge 
radius has reached roughly half of the radius Rq of the 
original droplets. We find that the volume of the liquid 
droplet stays constant during the whole coalescence pro- 
cess with a relative change of the order of 10~ 4 . We have 
chosen comparable material parameters (i.e., viscosity, 
surface tension and density contrast between inner and 
outer fluid) to the silicon oil used in Q, amounting to 
typical droplet radii of a few mm and time units between 
0.1 and 2 fis per l.u.. The viscous and inertial time scales 
are given by 



and 



tjRq 
a 



(12) 



(13) 



respectively, and range between 20 and 100 ms when ex- 
pressed in physical units. Detailed simulation parameters 
can be found in Table |U We cover a wide range of val- 
ues for surface tension, viscosity ratio (which is equal to 
the density ratio here), Ohnesorge and Reynolds number. 
The Ohnesorge number is given by 



Oh=^ = ^^ 



RoRe 



(14) 



and has been recently shown to be an important quantity 
separating different coalescence regimes [l7| . According 
to the phase-diagram proposed in [17| , the present values 
of the Ohnesorge number indicate that our simulations 
are located in the inertially-limited viscous regime. We 
will return to this point in sees. IIIIB and |V] 




As our LB simulations are based on a diffuse inter- 
face model [1, [25| , some care is needed in order to avoid 
spurious influences of the finite interface width on the 
determination of the bridge radius. Denoting by £ the 
characteristic length scale over which the liquid-vapor 
interface is diffuse in our simulations, we have to make 
sure that we only consider the regime where all our phys- 
ically interesting quantities are significantly larger than 
£, thereby approximating the sharp interface limit. This 
leads to the restrictions 5/£ ^> 1 for the bridge radius 
and rv/i; > 1 for the radius of curvature, r CT , of the neck 
(see Fig. [T]). While the first condition is easily fulfilled 
in our simulations, the latter one imposes a lower limit 
on the physically meaningful bridge radius. From the ge- 
ometric considerations of section IIV1 we have [eq. (|19l) ] 
r a = b 2 /2Rq, and thus arrive at the condition 



(15) 



With a typical value of £ = 4, we obtain b/Ro > 0.04 for 
R a = 5000, b/Ro > 0.15 for R = 400 and b/R > 0.3 for 
Ro = 100. 



III. RESULTS 
A. General aspects 

In Fig. a typical coalescence event in 2D as ob- 
served in our simulations is visualized. The upper inset 
of Fig. [5^ shows the evolution of the center-of-mass po- 
sition Ay cms = y cms (t) (0) of one of the droplets 
in dependence of the bridge radius b/Ro- For all our 
simulations, we find that the data is well described (in- 
cluding the proportionality constant) by a power-law 
Ay cms /Ro ~ 0.3(6/i?o) 3 ' 4 , represented by the solid line 
in the inset of Fig. 

We remark that if one plots Ay cms /Ro vs. t/r v rather 
than vs. b/Ro, a general power-law exponent of 1.7 is 
found, effectively in agreement with the predictions of 
the Stokesian theory of coalescence by Hopper In 
this representation, however, the prefactor turns out to 
be strongly dependent on the droplet size. The exponent 
3.4 that describes the dependence of Ay cms on b/Ro can 
be obtained from Hopper's theory by substituting the 
proper time-dependence of the bridge radius as observed 
in the present simulations (6 ~ t 1 / 2 , see below). We also 
note that, in contrast to predictions of 15] , we have never 
observed the presence of a vapor bubble around the neck. 

Fig. Hp shows the momentum field, pu, in the region 
close to the bridge typically observed in our simulations 
for a large density ratio (here, pl/pv — 100). We observe 
a strong mass transport from the interior of each drop 
towards the bridge, driven by the difference in Laplace 
pressure between the drop and negatively curved neck. 

Interestingly (see color field in Fig. [5p) , we find that 
the velocity field has a negative divergence (V • u < 0) in 
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TABLE I: Parameters for our simulations in three and two dimensions. The first column refers to the symbols used in the 
plots, r/L/rfv the liquid- vapor viscosity ratio (which is equal to the liquid- vapor density ratio), a the surface tension (in l.u.), 
r v and Ti are the viscous and inertial times (in l.u.), Oh is the Ohnesorge number [eq. (|14p ] and Re is the average Reynolds 
number (computed by averaging the maximum velocity in the system over the full time domain of the coalescence process). In 
the 2D case, the letters 'P and 's' below Oh refer to the simulations performed using a large (Ro = 5000) or small (Ro = 400) 
droplet radius. The values for t v and r,: are stated for the smaller radius. 
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FIG. 2: (a) Coalescence of two droplets initially at rest (in 2D). Contours represent the liquid- vapor interface and correspond to 
times 0, 0.01, 0.025, 0.05, 0.075, 0.15, 0.25, 0.375, 0.5 (in units of t v , from inner to outer). The bottom inset is a magnification of 
the bridge region. The upper inset shows the center-of-mass position, Ay cms = i/ C ms (t) — y C ms (0) of one of the droplets (measured 
along the symmetry axis) in a double-logarithmic representation. The solid line represents the prediction of Hopper [lj • (b) 
Momentum field pu at t ~ 0.075tj, in the region of the bridge for a typical simulation with a density ratio of Pl/pv = 100. 
The color-field represents the divergence of the velocity field, V • u (the darker the color the more negative the value; positive 
values are very small and not resolved by the color scaling, corresponding to white in the plot). 



the neck region, indicative of condensation (cf. 26] ) . The 
possibility of phase transition is in fact a characteristic 
feature of the present simulation model. The occurrence 
of condensation can be understood as a consequence of 
the Kelvin- (Gibbs- Thomson-) effect [27j and the fact that 
we study coalescence of liquid droplets that are initially 
in equilibrium with their vapor phase: while the droplets 
are separated, the interface curvature is positive through- 
out, implying a slightly elevated chemical potential corre- 
sponding to an increased saturation vapor density (reck- 



oned with respect to a flat interface) [28|. In contrast, 
the strongly negatively curved neck established after con- 
tact entails a local depression of the chemical potential, 
which in turn drives a condensation flux from the outer 
vapor regions towards the neck. 



For a more detailed study of this issue, we compare in 
Fig. [3] the different transport mechanism that contribute 
to the motion of the interface and their time dependence. 



5 



5. x i(r 5 

4. x 10" 5 
3. x 10" 5 
2. x 10" 5 
1.x 10" 5 



(a) 






i>Vp 


pV-u 







0. 

0.00 0.05 



0.10 0.15 

t/t v 



0.20 0.25 



1. x 10" 4 
8.X10 -5 
6. x 10" 5 
4. x 10 -5 

2. x 10" 5 

0. 



u-Vp 



pV-u 



(b) 



0.00 0.02 0.04 0.06 0.08 0.10 0.12 

t/ty 



FIG. 3: Comparison of the different mechanisms of mass transport [see eq. (|16p ] contributing to the interface motion near the 
neck, for a density ratio of (a) Pl/pv = 10 and (b) Pl/pv = 1000. Shown are the contributions from the advective term u- Vp 
and the divergence term pV • u, computed as a weighted average over the interface near the neck. 



The continuity equation [eq. Q], 

dtp = — pV • u — u • Vp , (16) 

implies that the interface can move both due to 
advection, represented by u ■ Vp, and evapora- 
tion/condensation, represented by pV • u. To compare 
the contribution of both terms, we compute them as a 
weighted average over the interface in a narrow region 
near the neck, i.e., J dr|V p(r)|s(r) / J dr|Vp|, where s(r) 
is either pV • u or u • Vp. As seen from Fig. [3J advection 
represents the dominant mechanism of interface motion 
in our system in all cases. As expected, the relative con- 
tribution form condensation increases at early time as 
the density ratio is reduced, since the compressibility of 
the model is enhanced. However, in the interesting late- 
time regime, we find that condensation is suppressed by 
almost an order of magnitude compared to advection, 
indicating that our model can be treated as effectively 
incompressible in this regime. 



B. Time-evolution of the bridge: 2D 

We now turn to a detailed investigation of the time 
evolution of the bridge radius in our simulations. Fig. [4^, 
contains the raw data of the time evolution of the bridge 
radius obtained from five different simulations in two di- 
mensions (for simulation parameters see Table . From 
the values of the Reynolds number (see inset) it can al- 
ready be inferred that the coalescence process is located 
in the viscous regime. [29| This is also corroborated by the 
(approximate) data collapse obtained when rescaling the 
time axis with the viscous time t v [eq. (1121) ]. as seen in 
Fig. 2}}. In contrast, when expressing the time in units 
of the inertial time Tj, no collapse is obtained. As the 
double-logarithmic representation of Fig. 0t shows, the 
bridge radius approximately follows a power-law b(t) ~ 
{t/T v ) a with an exponent a varying between 0.50 and 
0.58 within the individual curves. Note that there are 
no fitting parameters involved in obtaining these plots. 



Interestingly, this behavior is consistent with the predic- 
tion obtained from a simple scaling argument based on 
the incompressible Stokes equations [see section HV] . 

For comparison, Fig. [5] shows the evolution of the 
bridge radius obtained for droplets of radius i?o = 400, 
using the same set of simulation parameters as in Fig. S) 
Due to the smaller droplet radii, the evolution spans only 
about one decade in time. Nevertheless, by rescaling 
the time by the viscous time scale (Fig. [5^,) and plot- 
ting the data in logarithmic representation (Fig. \5jp) we 
again conclude that the droplets coalesce in the viscosity- 
dominated regime, following a W 2 -law in time. 

In order to examine the possible effect of grid reso- 
lution on our results, we have systematically increased 
the spatial resolutions at a constant Ohnesorge number 
[eq. ([11])] . Since the interface thickness is specific to the 
numerical scheme and does not have a direct counterpart 
in the corresponding macroscopic hydrodynamic equa- 
tions, we have kept in these studies the number of in- 
terface grid points constant so that a higher resolution 
also corresponds to a sharper interface. Results on the 
time evolution of the bridge radius obtained from these 
simulations are shown in Fig. [6l In plotting the data, we 
neglect the short initial transient which does not obey 
eq. (fT5)) . As seen from Fig. [51 upon rescaling b by Rq and 
the time t by the viscous time t v , all curves settle on a 
master curve ~ i ' 5 . We thus conclude that, already for 
the smallest investigated droplet sizes of Ro = 400, the 
underlying grid resolution is sufficiently high in order to 
resolve the correct hydrodynamic behavior. 

We remark that experiments on viscous coalescence of 
two-dimensional liquid lenses reported a linear time 
evolution for the bridge radius (similar to the 3D case), 
rather than a f^-scaling law. Based on a recent study 
of Paulsen et al. [I?} , one would expect our simulations 
to reside in the inertially-limited viscous regime, which 
is empirically characterized by a linear time evolution. 
For comparison, an analytical theory for Stokesian coa- 
lescence by Hopper [l4| predicts that b oc (t/T v )°- 857 in 
the range 0.0042 < b/Ro < 0.21 (note that this is an ap- 
proximation to the full solution provided by Hopper [14j , 
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FIG. 4: Viscous coalescence in 2D for two droplets of initial radius Ro = 5000 l.u. (a) Raw data for the time evolution of 
the bridge radius b (in units of the droplet radius Ro) as obtained from our simulations. The inset shows the evolution of 
the Reynolds number for one case. In (b), the data is shown when the time is expressed in terms of the viscous time r v . For 
comparison, in the inset, time is rescaled by the inertial time r*. In (c), the data is plotted in a double- logarithmic representation 
using the same rescaling of the axes as in (b). The dashed, dot-dashed and dotted lines represent power laws ~ t 1 ^ 2 , ~ f - 857 
and ~ t, respectively. See Table U for a legend to the different symbols. 



valid during the whole coalescence process). However, it 
is clearly visible from Fig. 2}d and c that there appears 
no sufficiently extended region where either a linear or 
i 857 -law can describe the data well. We will return to a 
discussion on possible reasons for these discrepancies in 
section [V] 



C. Time-evolution of the bridge: 3D 

Turning to coalescence in the three dimensions, we an- 
alyze our data in an analogous manner as in the 2D case. 
In Fig. [7^, we show the bridge radius evolution obtained 
from five representative simulation runs (see Tab. U for 
simulation parameters). Although simulation parame- 
ters are comparable to 2D, the Reynolds number remains 
roughly one order of magnitude smaller than in two di- 
mensions. As seen in Fig.[7b, expressing the time in units 
of the viscous time, t v , produces a reasonable scaling col- 
lapse of all data points onto a master curve. In contrast, 
no collapse is found when rescaling the time with the 
inertial time Tj (inset to Fig. [7b). This confirms, inde- 



pendently from the value of the Reynolds number, the 
dominance of viscosity in the coalescence process. As 
the double-logarithmic representation of Fig. |7fc shows, 
the data follows a power law b(t) ~ (i/r,,) 1 / 2 with rea- 
sonable accuracy. The power-law exponent is found to 
be only approximately equal to 1/2, varying in the range 
0.45 — 0.6 between the individual curves. Again, no fit- 
ting parameters are involved in the presentation of the 
data. Note that, due to computational limitations, the 
initial radii of the two droplets had to be chosen signif- 
icantly smaller than in 2D and, consequently, the time 
axis in Fig. [7] spans only roughly one decade. This range 
is nevertheless comparable to experiments [H, [l6[ . 

Similarly to the two-dimensional case, even when ac- 
knowledging the slight variation in the power-law ex- 
ponent, our results stand in contrast to experiments 
[1> Hi EH) where typically a linear time-dependence is 
observed (see, however, [1JJ). It is also interesting to 
compare our results to the analytical theory of three- 
dimensional coalescence of Eggers et al. [15j . which pre- 
dicts the bridge radius (in the region rh < b < Rq) to 
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FIG. 5: Viscous coalescence in 2D for two droplets of initial radius Ro = 400 l.u. Time is expressed in terms of the viscous 
time t v , while the bridge radius b is expressed in units of Ro- For comparison, in the inset, the time is rescaled by the inertial 
time n. In (a), the data is shown in a linear and in (b) in a double-logarithmic representation. The dashed and dotted lines 
represent power laws ~ t 1 / 2 and ~ t, respectively. See Table [J for a legend to the different symbols. 
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FIG. 6: Effect of grid refinement on the scaling behavior of the bridge radius b in 2D. For a fixed interface thickness of £ = 4 
(l.u.), we vary the drop radius while keeping the Ohnesorge number constant. The droplet radii are chosen as Ro = 400 (□), 
1600 (A), 5000 (•), 10000 (■) (all in l.u.). This corresponds to increasingly sharper interfaces: £/R = 0.02 (□), 0.005 (a), 
0.0016 (•), 0.0008 (■). Time is expressed in terms of the viscous time scale t v . In all cases, the density ratio is Pl/pv = 1000, 
the Ohnesorge number is 0.27 and the interface width as well as the initial separation of the droplets are kept constant at 4 
lattice units. 



evolve according to the differential equation 



b(t/r v ) 1 / 



b{th 



(17) 



Here, Co and a are constants that, in the original work 
15], are predicted as Co = 1 and a — 3 for an inviscid 
outer fluid while a = 3/2 for finite viscosity ratios (both 
valid for b/R < 1). For b < 0.03i? , the solution of 
eq. (|17|) follows an approximate scaling law of the form 

m 



b{t/T v ) 

Rq 



(18) 



In Fig. [5^,, the typical behavior of the numerical solution 
of eq. dT7|) as well as the approximate law (fT5j) is shown. 



Several interesting observations can be made: First of 
all, over many decades, the bridge radius as predicted 
by the analytical theory evolves significantly slower than 
it would be for a linear time-dependence. [30j] In fact, in 
the considered range, the behavior of eqs. (JTTJ) and (fTg|) 
can be well approximated by an effective power-law ~ t n , 
with an exponent n that is close to the value obtained by 
Hopper in the two-dimensional case [14j, n « 0.86 (the 
value of n increases for still earlier decades in time). At 
later times, the evolution of the bridge slows down and 
a time window appears where b approximately evolves 
proportionally to t 5 . 

Given the latter observation, it is tempting to fit the 
numerical solution of eq. (|17l) to our data for the bridge 
radii (Fig. [Hp). We find that reasonable agreement can 
only be obtained if a and cq are treated as fit parame- 
ters. In the inset to Fig. [Bp it is seen that the obtained 
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FIG. 7: Viscous coalescence in 3D. (a) Raw data for the time evolution of the bridge radius b (in units of the droplet radius 
Ro) as obtained from our simulations. The inset shows the evolution of the Reynolds number for one case. In (b), the data is 
shown when the time is expressed in terms of the viscous time t v . For comparison, in the inset, time is rescaled by the inertial 
time Ti. In (c), the data is plotted in a double- logarithmic representation using the same scaling of the axes as in (b). The 
dashed and dotted lines represent power laws ~ t 1 ' 2 and ~ t, respectively. See Table [T] for a legend to the different symbols. 



best-fit values for the scaling exponent a decrease with 
the liquid- vapor viscosity ratio rj^/rfv, settling around 
a value of 2 for large viscosity ratios - a trend which 
is opposite to the theoretical predictions of 15]. The 
parameter cq is found to be of the order of unity for 
all except the smallest viscosity ratios - a result which 
roughly agrees with the theoretical expectation, co = 1. 
We remark that, due to the limited range in time, it 
can presently not be decided whether our observation of 
a £ x / 2 -law in the data for three-dimensional coalescence 
(Fig. [7k) is a manifestation of a finite-size effect asso- 
ciated with the slowing-down of the coalescence process 
or whether it actually represents a true deviation from 
the typical time evolution in the Stokesian theory of coa- 
lescence [lj|, similar to the two-dimensional case. Here, 
additional simulations over a much larger range of time 
will be necessary to clarify this issue. This is reserved for 
future work. 



IV. SCALING CONSIDERATIONS 



In sec. IIII Al we observed that the contribution of con- 
densation to the motion of the bridge is much reduced 
compared to advection and, in particular, has no signifi- 
cant influence on the time-evolution of the bridge radius. 
Indeed, we found the scaling law b ~ (t/ry) 1 / 2 to be quite 
robust over a wide region of the parameter space, with 
the exponent generally being in the range 0.5 ±0.1. Moti- 
vated by this findings, we present here a simple derivation 
of this scaling law based on the incompressible Stokes 
equations - that is, neglecting evaporation or condensa- 
tion from the outset. 

Fig. [9] shows the basic situation we consider in our 
derivation. Rq is the initial radius of each of the two 
droplets, b the radius of the connecting bridge, and r a 
the radius of curvature of the meniscus. Due to rotational 
symmetry, it suffices to consider the problem in a plane 
that contains the conjoining line of the centers of the two 
droplets. By geometry, Rq = (b + r a ) 2 + (R — r a ) 2 , and 
thus, neglecting 2r 2 , we obtain the radius of curvature 



9 




FIG. 8: (a) Theoretical time evolution of the bridge radius according to the numerical solution of eq. (|17p (thin solid curve) 
and its asymptotic approximation [eq. (|18[) . thick solid curve], for a typical set of droplet parameters. At early times, the 
analytical predictions can be well approximated by an effective power-law b ~ t 86 (dot-dashed line). For comparison, also a 
linear power-law (b ~ t, dotted) is indicated, (b) Numerical solution of eq. (|17|) (solid curves) fitted to the data for the bridge 
radius b obtained from simulations in 3D (symbols). The inset shows that, in contrast to theoretical expectation [T^|, the fit 
parameter a decreases with increasing liquid-vapor viscosity ratio. Error bars are of the order of the symbol size. 



2i?o 



TTnT ' (19) 



where the last approximation is justified since we focus 
on b <C Rq- It should be remarked that a similar re- 
sult is also used in [H, but a number of previous works 
[lil [HI found r a ~ b 3 . To obtain a relation for the time 
dependence of the bridge radius, we apply scaling argu- 
ments similar in spirit to [ g,ll5l Il8| ] to the incompressible 
Navier-Stokes equations LU 




p(d t + u • V)u = -Vp - 



r?V 2 u, 



(20) 



where p is the density, u the fluid velocity and p the 
local pressure. We shall neglect a possible rigid trans- 
lation of the coalescing droplets (although, in principle, 
this effect can be relevant [Tt]]) and also ignore the 



influence of the vapor on the dynamics. Focusing now 
on the axis along the direction of the bridge radius, the 
flow velocity u near the meniscus is, by continuity, deter- 
mined by the motion of the bridge alone, u = b. Assum- 
ing that the pressure varies between zero at the center of 
the bridge and the Laplace pressure pl — —a/r a at the 
curved meniscus, we can approximate Vp ~ pl /b. Fi- 
nally, we shall assume that the velocity varies smoothly 
(e.g., parabolically) between the center of the bridge 
(where u ss 0) and the interface, allowing us to estimate 
V 2 m ~ -u/b 2 . 

In the viscous regime, eq. l|2"0)) becomes rj\7 2 u = Vp^ 
in the steady-state and above simplifications lead to 



2R <t 1 
r] b 



(21) 



The solution of this differential equation is easily written 
down as 




FIG. 9: Geometry of the bridge region of two coalescing 
droplets as considered in our scaling theory. Ro is the droplet 
radius, b the radius of the connecting bridge and r a the radius 
of curvature of the meniscus. 



b(tf-bl = Rl^>, (22) 

T 

where to is a suitably chosen initial time, bg = b(to) 
and the characteristic time scale r is (up to a numeri- 
cal prefactor) given by the viscous time t„/4 = tjRq/Ao- 
[cq. (fl~2"j)]. For completeness, we state also the result for 
the inertial regime, where the advection term p(u • V)u 
dominates over the viscous stress. Here, the same so- 
lution as in eq. (|2"2"|) is obtained, but with r replaced by 
the inertial time Tj/v8 = \J ' pR^/Sa [eq. ([X5|], It must be 
emphasized here that the above results hold independent 
of the spatial dimension. 

Neglecting the integration constants to and b(to) in 
(|22p. which is justified except in the early stages of the 
evolution, the growth-law for the bridge radius can be 
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expressed as 



Ro 



1/2 



(23) 



In the inertial regime, the above relation is well known 
@ and has been confirmed by experiments 0, H, 0, H, EH ■ 
Remarkably, our simple scaling theory yields a power-law 
with the same exponent 1/2 also in the viscous regime, 
albeit with a different characteristic time scale. We note 
that relation (|23l) agrees with Frenkel's classic result on 
viscous sintering [12j, where, however, it was obtained 
under the assumption that the two coalescing droplets 
reduce their surface free energy by moving closer together 
as a whole, without considering explicitly the dynamics 
of the bridge. 

Of course, the scaling arguments leading to eq. (|2"3")l 
are not rigorous and therefore have to be taken with care. 
For instance, approximating the Laplacian of the velocity 
at the bridge as V 2 it ~ u/b n with an arbitrary power n 
allows one to write down a generalization of eq. ([2~Tj) as 



n — 3 



6 = ^fA 



Obviously, the choice n — 2 corresponds to eq. (f2"B1 , sup- 
ported by our simulations, while n = 4 would lead to a 
linear time dependence of b, in agreement with experi- 
ments. 

In sec. IHII we have shown that the time evolution of the 
bridge radius observed in our simulations is consistent 
with a scaling relation of the form of eq. (f23|) . without 
invoking any fit parameters. Assuming for the moment 
the validity of the present derivation, we can make use of 
the freedom to choose the initial value of the time to as 
allowed by the solution (|2"2"j) of our scaling ansatz. The 
bridge radius data transformed according to the relation 
(22), where bo = b(to) is fixed by the choice of to, is shown 
in Fig. [TOa.b in the two-dimensional and in Fig. [TUb in 
the three-dimensional case. Consistent with the results 
in sec. IIII| we find reasonably good scaling collapse of all 
data points if time axis is scaled by the viscous time, and 
no collapse if the inertial time is used (see insets). 



V. SUMMARY 

In this work, we have investigated the coalescence of 
two identical resting droplets in the viscous regime in two 
and three dimensions. From Lattice Boltzmann simula- 
tions of the full isothermal Navicr- Stokes equations for 
a non-ideal fluid, we find that the time evolution of the 



bridge radius can be well described by a i^-scaling law - 
in striking similarity to what is generally found to hold in 
the inertial regime, albeit with a different characteristic 
time constant. The velocity field shows a non-zero diver- 
gence at the neck of the coalescing droplets, indicating 
the presence of condensation. However, over the investi- 
gated parameter range, condensation is found to give a 
negligible contribution to the motion of the interface at 
the neck, the dominant process being advection. Guided 
by this observation, we have shown that our simulation 
results can be rationalized through simple scaling argu- 
ments applied to the incompressible Stokes equations [see 
eq. ((22])]. 

Our findings differ markedly from recent experiments 
0) IS B EH ! which report a linear time-evolution of the 
bridge radius in the viscous regime both in two and three 
dimensions. In two dimensions, our findings are also 
clearly different from analytical theories [14|, which pre- 
dict a power-law b ~ t 86 . However, the center-of-mass 
motion of each droplet is found to be described by a 
power-law ~ t 1 ' 7 , with an exponent that, surprisingly, 
agrees well with the theoretical prediction [3]. Due to 
the limited range in time, no definite conclusions con- 
cerning possible agreement with the theory for three- 
dimensional Stokesian coalescence [l5| can be drawn at 
present. It is also useful to note that, in the framework 
of the phase-diagram of [171 ] , our simulations would not 
be located in the true Stokes regime, which applies to 
droplets of much larger viscosity than presently used, but 
instead in the inertially-limited viscous regime. However, 
this does not resolve the above mentioned discrepancies 
since most experiments also reside in this region [j, d, 0] ■ 
Furthermore, the time-evolution of the bridge radius in 
the inertially-limited viscous regime has only empirically 
been found to follow a linear time evolution, with a the- 
oretical derivation of this result lacking so far. Presently, 
we have no explanation for the discrepancies between our 
work and existing theories or experiments. This is a sub- 
ject of future work. 
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